9.10 Integration of Bayesian Inference Techniques with Mathematical Modeling

The credibility of the scientific methodology of numerical models and their adequacy to form the basis of public policy decisions have been frequently challenged. The first part of this chapter aims to address the issue of model reliability by evaluating the current state of aquatic biogeochemical modeling. We provide evidence that there is still considerable controversy among modelers and the resource managers about how to develop, evaluate, and interpret mathematical models. Our arguments are that (1) models are not always developed in a consistent manner, clearly stated purpose, and predeter­ mined acceptable model performance level, and (2) the potential users select models without properly assessing their technical value. The second part of this presentation argues that the development of novel methods for rigorously assessing the uncertainty underlying model predictions should be a top priority of the modeling community. Striving for novel uncertainty analysis tools, we introduce Bayesian calibration of process-based models as a methodological advancement that warrants consideration in aquatic ecosystem research. This modeling framework combines the advantageous features of both process-based and statistical approaches, that is, mechanistic understanding that remains within the bounds of data-based parameter estimation. The incorporation of mechanism improves the confidence in predictions made for a variety of conditions, whereas the statistical methods provide an empirical basis for parameter value selection and allow for realistic estimates of predictive uncertainty. Other advantages of the Bayesian approach include the ability to sequentially update beliefs as new knowledge is available, and the consistency with the scientific process of progressive learning and the policy practice of adaptive management. Finally, we illustrate some of the anticipated benefits from the Bayesian calibration framework, well suited for stakeholders and policy makers when making environmental management decisions, using the Hamilton Harbour – a eutrophic system in Ontario, Canada – as a case study. 9.10.1 Evaluation of the Current State of Aquatic Biogeochemical Modeling: Where Are We? Mechanistic aquatic biogeochemical models have formed the scientific basis for environmental management decisions by providing a predictive link between management actions and ecosystem response. They have provided an important tool for elucidating the interactions between climate variability and the Treatise on Estuarine and Coastal Science, 2011, Vol.9 carbon cycling in the oceans, and thus for assessing the pace and impacts of climate change (Doney, 1999; Franks, 2002). Acknowledging their central role in aquatic ecosystem research, several compelling questions arise, such as: What is the capacity of the current models to simulate the dynamics of coastal and estuarine ecosystems? How carefully do modelers develop their models? How rigorously do we assess what a model can or cannot predict? Arhonditsis and Brett (2004) attempted to 173 , 173-192, DOI: 10.1016/B978-0-12-374711-2.00910-4 Author's personal copy 174 Integration of Bayesian Inference Techniques with Mathematical Modeling answer some of these questions by reviewing 153 modeling studies published in the literature between 1990 and 2002. Their hypothesis was that the sizable number of aquatic eco­ system modeling studies, which successfully passed the scrutiny of the peer-review process along with the experience gained from addressing an extent of management problems, can objectively reveal the systematic biases, methodological inconsistencies, and common misconceptions characterizing the field of aquatic ecosystem modeling. Indeed, despite the heterogeneity of the modeling studies examined with respect to model complexity, type of ecosystem modeled, spatial and temporal scales, and model development objectives, this study was able to detect statistically significant trends of the model performance and to pinpoint methodological omis­ sions in the current modeling practice. The first interesting finding was the absence of systematic goodness-of-fit assessment of the original models, that is, plots in which simulated values were visually compared with observed data were only presented for 16.8% of the model endpoints, and even fewer (1.3%) were the cases in which thorough statistical examination of the error was reported. In the cases in which measures of fit or comparison plots were presented, Arhonditsis and Brett (2004) independently assessed state variable performance as expressed by the relative error (RE = Σ |observed values – simulated values|/Σ observed values) and the coefficient of determination (r) (Table 1). It was found that temperature and dissolved oxygen had the lowest RE (median < 10%) and the highest r values (the respective medians were 0.93 and 0.70). The typical limiting nutrient forms (NO3, NH4, PO4, and Si) in freshwater and oceanic ecosystems along with the phytoplankton biomass had intermediate fit, with median r values varying from 0.40 to 0.60 and the median RE lying around the 40% level. Zooplankton dynamics were characterized by the highest RE (70%) and the widest range of r (interquartile range ∼0.8) and Table 1 Performance of the aquatic biogeochemical models for the study Percentile Temperature Dissolved oxygen Nitrate Ammonium 10th r 2 0.42 0.34 0.10 0.05 RE(%) 2 4 8 18 20th r 2 0.62 0.52 0.37 0.13 RE(%) 4 7 18 30 30th r 2 0.81 0.58 0.47 0.18 RE(%) 5 8 26 34 40th r 2 0.92 0.62 0.56 0.29 RE(%) 5 10 32 40 50th r 2 0.93 0.70 0.68 0.39 RE(%) 7 12 36 48 60th r 2 0.95 0.78 0.79 0.44 RE(%) 7 14 44 55 70th r 2 0.96 0.86 0.84 0.57 RE(%) 9 17 57 65 80th r 2 0.97 0.88 0.91 0.78 RE(%) 90th r 2 11 0.98 19 0.92 68 0.95 77 0.89 RE(%) 15 22 88 101 100th r 2 0.99 1.00 1.00 0.99 RE(%) 25 31 554 206 Coefficient of determination (r ) and relative error (RE%) values for: temperature, dissolved ox Adapted from Arhonditsis, G.B., Brett, M.T., 2004. Evaluation of the current state of mechanis Treatise on Estuarine and Coastal Science, 2011, Vol.9, 1 RE (interquartile range ∼85%) values. Similarly, bacteria were also poorly modeled (median r value <0.06), indicating that the performance of existing mechanistic aquatic biogeochem­ ical models declines as we move from physical–chemical to biological components of planktonic systems. On a positive note, it was found that these results were obtained without the introduction of a major ‘calibration bias’, that is, in the process of maximizing the fit for a specific state variable (usually phytoplankton biomass), the modelers did not seem to com­ promise on the fit for other state variables (such as limiting nutrient concentrations or herbivorous zooplankton biomass). Arhonditsis and Brett (2004) also assessed the effects of model complexity (expressed as the number of state variables), spatial dimension (from zeroto three-dimensional models), simulation period (from days to decades), and ecosystem type on model performance. The study reported a positive correla­ tion between the number of state variables and the RE values for the different model outputs (r = 0.219, p < 0.001). This (counterintuitive) positive trend was even stronger when con­ sidering the RE values for phytoplankton (r = 0.248, p = 0.003) and zooplankton (r = 0.626, p < 0.001) biomass, suggesting that more complex models usually result in slightly poorer model performance. It should be noted, however, that the majority of the complex models considered in this analysis belonged to the European Regional Seas Ecosystem Model (ERSEM) family, and therefore the reported complexity– performance relationship was influenced by the development purposes, modeled environments and practices, followed by this particular family of models (Baretta et al., 1995). Similarly, a (very weak) positive correlation was found between the dura­ tion of the simulation period and the state variable RE values (r = 0.098, p = 0.022), indicating that longer simulations are also increasing model misfit. Marginally significant correla­ tions also exist between the spatial complexity of the models and their (RE values) performance trends (r = 0.104, p = 0.015).